classdef Mesh < handle
    %UNTITLED2 Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        numel
        nnod
        ndm
        nen
        nodes;
        eles;
        dim
        %para que sea mas facil recuperar los datos para la ventana
        Lx
        Nx
        Ly
        Ny
        refx
        refy
        type
        
        K_diffusion %matriz de rigidez (incluye cb robin si la hay)
        b_diffusion %flujo difusivo (incluye cb robin si la hay, y flujo si lo hay)
        K_convection %matriz convectiva
        b_source %fuente
        K_reaction %matriz reactiva
        K_temporal %matriz de masa temporal
        b_temporal %carga temporal
        
        K_global %matriz global
        b_global %vector derecho global
    end
    
    methods
        function [obj] = Mesh(obj)
            
        end
        
        %una vez elegidos los parametros a utilizar y creadas las funciones
        %custom, se procede a empaquetar los datos para resolver el
        %problema.
        function [obj] = setProblem(obj)
            global data
            param = data.param;
            t = param.ttime;
            addpath(param.direc);
            xx = get(obj.nodes,'x');
            
            %verifico si el coeficiente es constante o custom
            if get(param,'cust_diff')
                Kvalues = custom_diffusion(xx,t);
                set(data.mesh.nodes,'k',num2cell(Kvalues));
            else
                set(data.mesh.nodes,'k',{get(param,'diffusion')});
            end
            
        end
        
        function [obj] = set(obj,Lx,Nx,refx,Ly,Ny,refy,type)
            if nargin == 4 
                % caso 1d - estructurado
                obj.Lx=Lx;
                obj.Nx=Nx;
                obj.refx=refx;
                obj.dim=1;
                obj.numel = Nx;
                obj.nnod = Nx+1;
                obj.ndm = 1;
                obj.nodes(1,Nx+1) = Node;
                obj.eles(1,Nx) = Element;
                
                xx = refh(refx,Nx+1)*Lx;
                
                %haciendo un for tradicional (N=10000 4.59s)
                %tic
                %for i=0:Nx
                %    set(obj.nodes(i+1),'x',i*Lx/Nx);
                %end
                %toc
                
                %vectorizado (N=10000 0.02s)
                set(obj.nodes,{'x'},num2cell(xx)); %necesito convertirlo a cellarray
                if Nx==2
                    set(obj.eles,{'icone'},num2cell([ [1:Nx]' [2:Nx+1]'],[2]));
                else
                    set(obj.eles,{'icone'},num2cell([ [1:Nx]' [2:Nx+1]'],[Nx,2])); %necesito convertirlo a cellarray
                end
                
                figure(1);clf;plot(cell2mat(get(obj.nodes,'x')),zeros(obj.nnod,1),'*');
                
            elseif nargin == 8
                obj.Lx=Lx;
                obj.Nx=Nx;
                obj.refx=refx;
                obj.Ly=Ly;
                obj.Ny=Ny;
                obj.refy=refy;
                obj.type=type;
                obj.dim=2;
                
                [xnod,icone] = mesh_square(Lx,Nx,Ly,Ny,refx,refy,type);
                
                [obj.nnod,obj.ndm] = size(xnod);
                [obj.numel,obj.nen] = size(icone);
                
                obj.nodes(1,obj.nnod) = Node;
                obj.eles(1,obj.numel) = Element;
                
                %vectorizado (N=10000 0.02s)
                set(obj.nodes,{'x'},num2cell(xnod,[obj.nnod,2])); %necesito convertirlo a cellarray
                set(obj.eles,{'icone'},num2cell(icone,[obj.numel,2])); %necesito convertirlo a cellarray
                
                figure(1);clf;pltmsh(cell2mat(get(obj.nodes,'x')),cell2mat(get(obj.eles,'icone')));axis equal;
            end
        end

        function [obj] = assemble(obj)
            global data;
            obj.K_global = sparse([],[],[],obj.nnod,obj.nnod);
            obj.b_global = zeros(obj.nnod,1);
            %si es custom o el valor constante es distinto cero -> se resuelve
            if get(data.param,'cust_diff') || get(data.param,'diffusion')
                obj.assemble_diffusion();
                obj.K_global = obj.K_global + obj.K_diffusion;
                obj.b_global = obj.b_global + obj.b_diffusion;
            end
            
            if get(data.param,'cust_conv') || get(data.param,'convection_x')  || get(data.param,'convection_y')
                obj.assemble_convection();
                obj.K_global = obj.K_global + obj.K_convection;
            end
            
            if get(data.param,'cust_sour') || get(data.param,'source')
                obj.assemble_source();
                obj.b_global = obj.b_global + obj.b_source;
            end
            
            if get(data.param,'cust_reac') || get(data.param,'reaction')
                obj.assemble_reaction();
                obj.K_global = obj.K_global + obj.K_reaction;
            end
            
            npatchs = length(data.patchs);
            
            for ipatch=1:npatchs
                if(strcmp(data.patchs(ipatch).type,'DIRICHLET'))
                    nn = cell2mat(data.patchs(ipatch).nodes);
                    Tref = cell2mat(data.patchs(ipatch).Tref);
                    obj.K_global(nn,:) = 0;
                    dd = diag(obj.K_global);
                    dd(nn) = 1;
                    obj.K_global = obj.K_global - diag(diag(obj.K_global)) + diag(dd);
                    obj.b_global(nn) = Tref*ones(size(nn,1),1);
                end
            end
            
            
        end
        
        function [obj] = assemble_diffusion(obj)
            
                row = [];
                col = [];
                coe = [];
                obj.b_diffusion = zeros(obj.nnod,1);
                %loop sobre cada elemento
                for iele = 1:obj.numel
                    obj.eles(iele).assemble_diffusion();
                    icone = get(obj.eles(iele),'icone');
                    K_ele = get(obj.eles(iele),'K_diffusion');
                    b_ele = get(obj.eles(iele),'b_diffusion');
                    [nen] = size(icone,2);
                    for inod=1:nen
                        for jnod=1:nen
                            row = [row;icone(inod)];
                            col = [col;icone(jnod)];
                            coe = [coe;K_ele(inod,jnod)];
                        end
                        obj.b_diffusion(icone(inod)) = obj.b_diffusion(icone(inod)) + b_ele(inod);
                    end
                end
                obj.K_diffusion = sparse(row,col,coe);                             
                
        end
        
        
        function [obj] = assemble_convection(obj)
        end
        
        function [obj] = assemble_reaction(obj)
            row = [];
                col = [];
                coe = [];
                %loop sobre cada elemento
                for iele = 1:obj.numel
                    obj.eles(iele).assemble_reaction();
                    icone = get(obj.eles(iele),'icone');
                    K_ele = get(obj.eles(iele),'K_reaction');
                    [nen] = size(icone,2);
                    for inod=1:nen
                        for jnod=1:nen
                            row = [row;icone(inod)];
                            col = [col;icone(jnod)];
                            coe = [coe;K_ele(inod,jnod)];
                        end
                    end
                end
                obj.K_reaction = sparse(row,col,coe);
        end
        
        function [obj] = assemble_source(obj)
                obj.b_source = zeros(obj.nnod,1);
                %loop sobre cada elemento
                for iele = 1:obj.numel
                    obj.eles(iele).assemble_source();
                    icone = get(obj.eles(iele),'icone');
                    b_ele = get(obj.eles(iele),'b_source');
                    [nen] = size(icone,2);
                    for inod=1:nen
                        obj.b_source(icone(inod)) = obj.b_source(icone(inod)) + b_ele(inod);
                    end
                end
        end
        
        function [obj] = assemble_temporal(obj)
        end
        
        function [obj] = solve(obj)
            phi = obj.K_global\obj.b_global;
            set(obj.nodes,{'state'},num2cell(phi));  
            xx=cell2mat(get(obj.nodes,{'x'}));
            if data.param.ndm==1
                figure(2);clf;plot(xx,phi,'b');
            else
                icone=cell2mat(get(obj.eles,{'icone'}));
                figure(2);clf;view2d(xx,icone,phi);shading interp; colorbar;
            end
            %xx=[0:0.001:1];yy=xx.^2/2+0.5*xx;hold on;plot(xx,yy,'r');legend('FEM','analitica')
        end
        
    end
        
            
end

